*========================================================================
*========================================================================
*======== BEGIN OUTPUT.do ===============================================
*========================================================================
*========================================================================
* Reproduction of figures and tables in
* "Who Values Future Energy Savings? Evidence from American Drivers"
* by Arik Levinson & Lutz Sager
clear all
use "VEH1.dta", clear

*========================================================================
**# Figure 1
*========================================================================
* Annual fuel cost differences, gas versus hybrid
	local XTITLE "Annual Cost Difference  ($, real 2017)"
	preserve
	keep if AnnualSavings_real < 1500

	twoway (histogram AnnualSavings_real if HybridVehicle==1,  color(olive_teal) lcolor(teal)  width(35) percent ) ///
	       (histogram AnnualSavings_real if HybridVehicle==0,  fcolor(none) lcolor(black)  width(35) percent ) ///
		   (kdensity AnnualSavings_real if HybridVehicle==1,  lcolor(green) yaxis(2) yscale(off range(0 0.002) axis(2)) ) ///
	       (kdensity AnnualSavings_real if HybridVehicle==0,  lcolor(black) yaxis(2) yscale(off range(0 0.002) axis(2)) ) ///
		    , ///
		   graphregion(color(white)) /*title("Hybrid and gas-powered vehicles (all hybrid/gas pairs)")*/ ///
			xtitle("`XTITLE'")  ytitle("Density (%)") ///
			legend( label(1 "Hybrid") label(2 "Gas-powered")  ///
			     position(2) bmargin(medium) symxsize(medium) size(small) rows(1) ring(0) region(lwidth(none)) ///
			     order(2 1 )) 	
				 
 		   graph export "output/Figure_01.png", as(png) width(1000) replace	
	restore


*====================================
**# Figure 2
*====================================
* Stylized figure not using any data.

*========================================================================
**# Figure 3
*========================================================================
	local XTITLE "Annual Cost Difference  ($, real 2017)"
	preserve
		keep if AnnualSavings_real < 1500
		twoway histogram AnnualSavings_real if HybridVehicle==0,  fcolor(none) lcolor(black)  width(35) frequency yscale(range(0 2000)) yla(0(500)2000) xline(300) ///
		 text(500 150 "A", place(e) color(red) box size(vlarge)) ///
		 text(200 500 "B", place(e) color(red) box size(vlarge)) ///
		   graphregion(color(white))  ///
			xtitle("`XTITLE'")  ytitle("Frequency") ///
				legend( label(1 "Gas-powered") ///
			     position(2) bmargin(medium) symxsize(medium) size(small) rows(1) ring(0) region(lwidth(none)) )	 
 		   graph export "output/Figure_03a.png", as(png) width(1000) replace	
	restore	
	local XTITLE "Annual Cost Difference  ($, real 2017)"
	preserve
		keep if AnnualSavings_real < 1500
		twoway histogram AnnualSavings_real if HybridVehicle==1,  color(olive_teal%50) lcolor(teal)  width(70) frequency yscale(range(0 2000)) yla(0(500)2000) xline(300) ///
		 text(120 150 "C", place(e) color(red) box size(vlarge)) ///
		 text(120 400 "D", place(e) color(red) box size(vlarge)) ///
		   graphregion(color(white))  ///
			xtitle("`XTITLE'")  ytitle("Frequency") ///
				legend( label(1 "Hybrid") ///
			     position(2) bmargin(medium) symxsize(medium) size(small) rows(1) ring(0) region(lwidth(none)) )	 
 		   graph export "output/Figure_03b.png", as(png) width(1000) replace	
	restore

		
*========================================================================
**# Table 2
*========================================================================
* Share of HAS and SHOULD (EPA 14yrs + 3% / 7% cutoffs, 2009/2017 NHTS)
		* Ratio of fuel savings based on fuel price in purchase year and msrp
		 gen FUELSAVING_vs_PREMIUM = FUELSAVING_vs_PREMIUM_initial		 
		*---------- IDENTIFY WHO SHOULD / SHOULD NOT ----------------		
		* APPROACH1_: cutoffs are calculated for an average 14 years of discounted fuel savings (at either 3% and 7% discounting)
			scalar APPROACH1_HI_cutoff = (1/8.75)		/* 14 years discounted at 7% */
			scalar APPROACH1_LO_cutoff = (1/11.30)		/* 14 years discounted at 3% */
			// identify those drivers whose cost-minimizing choice would have been the hybrid (SHOULD)
				gen APPROACH1_HI_SHOULD = .
				replace APPROACH1_HI_SHOULD = 0
				replace APPROACH1_HI_SHOULD = 1 if FUELSAVING_vs_PREMIUM > APPROACH1_HI_cutoff
				replace APPROACH1_HI_SHOULD = . if FUELSAVING_vs_PREMIUM ==.
				gen APPROACH1_LO_SHOULD = .
				replace APPROACH1_LO_SHOULD = 0
				replace APPROACH1_LO_SHOULD = 1 if FUELSAVING_vs_PREMIUM > APPROACH1_LO_cutoff
				replace APPROACH1_LO_SHOULD = . if FUELSAVING_vs_PREMIUM ==.
			// identify those drivers whose actual choice (HAS) differs from the cost-minimizing choice (SHOULD)
				gen APPROACH1_HI_HYBRID_WRONG = .
				replace APPROACH1_HI_HYBRID_WRONG = 0 if HAS == 1 & APPROACH1_LO_SHOULD !=.
				replace APPROACH1_HI_HYBRID_WRONG = 1 if HAS == 1 & APPROACH1_HI_SHOULD == 0
				gen APPROACH1_HI_NONHYBRID_WRONG = .
				replace APPROACH1_HI_NONHYBRID_WRONG = 0 if HAS == 0 & APPROACH1_LO_SHOULD !=.
				replace APPROACH1_HI_NONHYBRID_WRONG = 1 if HAS == 0 & APPROACH1_HI_SHOULD == 1
				gen APPROACH1_LO_HYBRID_WRONG = .
				replace APPROACH1_LO_HYBRID_WRONG = 0 if HAS == 1 & APPROACH1_LO_SHOULD !=.
				replace APPROACH1_LO_HYBRID_WRONG = 1 if HAS == 1 & APPROACH1_LO_SHOULD == 0
				gen APPROACH1_LO_NONHYBRID_WRONG = .
				replace APPROACH1_LO_NONHYBRID_WRONG = 0 if HAS == 0 & APPROACH1_LO_SHOULD !=.
				replace APPROACH1_LO_NONHYBRID_WRONG = 1 if HAS == 0 & APPROACH1_LO_SHOULD == 1
				capture drop SHOULD
				gen SHOULD = APPROACH1_LO_SHOULD
			
		preserve
			keep if SHOULD !=.
			outreg2 HAS APPROACH1_LO_SHOULD using "output/Table_02.xls", replace cross noaster
		restore
		preserve
			keep if SHOULD !=.
			outreg2 HAS APPROACH1_HI_SHOULD using "output/Table_02.xls", append cross noaster
		restore		

	
		
*========================================================================
**# Figure 4
*========================================================================
	
				gen HAS_HYBRID_BUT_BETTER_GAS = 0 if HAS == 1 & SHOULD == 1
				replace HAS_HYBRID_BUT_BETTER_GAS = 1 if HAS ==1 & SHOULD == 0
				gen HAS_GAS_BUT_BETTER_HYBRID = 0 if HAS ==0  & SHOULD == 0
				replace HAS_GAS_BUT_BETTER_HYBRID = 1 if HAS ==0 & SHOULD == 1
				sum HAS_HYBRID_BUT_BETTER_GAS HAS_GAS_BUT_BETTER_HYBRID

				gen SHOULD_HYBRID_BUT_HAS_GAS = 0 if SHOULD == 1
				replace SHOULD_HYBRID_BUT_HAS_GAS = 1 if HAS == 0 & SHOULD == 1
				gen SHOULD_GAS_BUT_HAS_HYBRID = 0 if SHOULD == 0
				replace SHOULD_GAS_BUT_HAS_HYBRID = 1 if HAS == 1 & SHOULD == 0
				sum SHOULD_HYBRID_BUT_HAS_GAS SHOULD_GAS_BUT_HAS_HYBRID
 
		* Panel (a): Shares by Income group. NHTS 2017 only.
				preserve
					keep if NHTSwave==2017
					collapse (mean) SHOULD_GAS_BUT_HAS_HYBRID SHOULD_HYBRID_BUT_HAS_GAS (semean) SD_SHOULD_GAS_BUT_HAS_HYBRID=SHOULD_GAS_BUT_HAS_HYBRID SD_SHOULD_HYBRID_BUT_HAS_GAS=SHOULD_HYBRID_BUT_HAS_GAS, by(hhfaminc_2017_groups)
					replace SHOULD_GAS_BUT_HAS_HYBRID = (SHOULD_GAS_BUT_HAS_HYBRID ) * 100
					replace SHOULD_HYBRID_BUT_HAS_GAS = (SHOULD_HYBRID_BUT_HAS_GAS ) * 100
					replace SD_SHOULD_GAS_BUT_HAS_HYBRID = (SD_SHOULD_GAS_BUT_HAS_HYBRID ) * 100
					replace SD_SHOULD_HYBRID_BUT_HAS_GAS = (SD_SHOULD_HYBRID_BUT_HAS_GAS ) * 100
					gen hi1=SHOULD_GAS_BUT_HAS_HYBRID+1.96*SD_SHOULD_GAS_BUT_HAS_HYBRID
					gen hi2=SHOULD_HYBRID_BUT_HAS_GAS+1.96*SD_SHOULD_HYBRID_BUT_HAS_GAS
					gen lo1=SHOULD_GAS_BUT_HAS_HYBRID-1.96*SD_SHOULD_GAS_BUT_HAS_HYBRID
					gen lo2=SHOULD_HYBRID_BUT_HAS_GAS-1.96*SD_SHOULD_HYBRID_BUT_HAS_GAS
					label variable SHOULD_GAS_BUT_HAS_HYBRID "Should Gas, but has Hybrid (actual)"
					label variable SHOULD_HYBRID_BUT_HAS_GAS "Should Hybrid, but has Gas (actual)"
					twoway (connected SHOULD_HYBRID_BUT_HAS_GAS hhfaminc_2017_groups, lwidth(1) lcolor(red) mcolor(red) msymbol(Oh) msize(large) lwidth(1) xlabel(1 2 3 4 5 6, valuelabel labsize(medsmall)) ylabel(,labsize(large)))  ///
							(connected SHOULD_GAS_BUT_HAS_HYBRID hhfaminc_2017_groups, lpattern(-) lcolor(green) mcolor(green) msymbol(D) msize(large) lwidth(1) xlabel(1 2 3 4 5 6, valuelabel labsize(medsmall)) ylabel(,labsize(large))) ///
							(rarea hi1 lo1 hhfaminc_2017_groups, lcolor(grey%30) color(grey%30)) ///
							(rarea hi2 lo2 hhfaminc_2017_groups, lcolor(grey%30) color(grey%30)) ///
							, ///
					   graphregion(color(white)) /*title("Likely Mistakes by Demographic")*/ ///
						xtitle("Income group", size(vlarge))  ytitle("Share of group (in %)", size(vlarge)) ///
						legend( label(1 "p(Gas | High Expenses)") label(2 "p(Hybrid | Low Expenses)")  ///
							 position(3) bmargin(medium) symxsize(medium) size(medium) rows(1) ring(0) region(lwidth(none)) ///
							 order(2 1 )) 	
					graph export "output/Figure_04a.png", as(png) width(1000) replace				
			restore

		* Panel (b): Shares by Income group. NHTS 2017 only. [DE-MEANED]
			preserve
				quietly capture gen HASGAS = 1-HAS
				quietly sum HAS
				scalar mean = r(mean)
				replace SHOULD_GAS_BUT_HAS_HYBRID = (SHOULD_GAS_BUT_HAS_HYBRID - mean) * 100
				quietly sum HASGAS
				scalar mean = r(mean)
				replace SHOULD_HYBRID_BUT_HAS_GAS = (SHOULD_HYBRID_BUT_HAS_GAS - mean) * 100
				keep if NHTSwave==2017
				collapse (mean) SHOULD_GAS_BUT_HAS_HYBRID SHOULD_HYBRID_BUT_HAS_GAS (semean) SD_SHOULD_GAS_BUT_HAS_HYBRID=SHOULD_GAS_BUT_HAS_HYBRID SD_SHOULD_HYBRID_BUT_HAS_GAS=SHOULD_HYBRID_BUT_HAS_GAS, by(hhfaminc_2017_groups)
				gen hi1=SHOULD_GAS_BUT_HAS_HYBRID+1.96*SD_SHOULD_GAS_BUT_HAS_HYBRID
				gen hi2=SHOULD_HYBRID_BUT_HAS_GAS+1.96*SD_SHOULD_HYBRID_BUT_HAS_GAS
				gen lo1=SHOULD_GAS_BUT_HAS_HYBRID-1.96*SD_SHOULD_GAS_BUT_HAS_HYBRID
				gen lo2=SHOULD_HYBRID_BUT_HAS_GAS-1.96*SD_SHOULD_HYBRID_BUT_HAS_GAS
				label variable SHOULD_GAS_BUT_HAS_HYBRID "Should Gas, but has Hybrid (pp. difference)"
				label variable SHOULD_HYBRID_BUT_HAS_GAS "Should Hybrid, but has Gas (pp. difference)"
				twoway (connected SHOULD_HYBRID_BUT_HAS_GAS hhfaminc_2017_groups, lwidth(1) lcolor(red) mcolor(red) msymbol(Oh) msize(large) lwidth(1) xlabel(1 2 3 4 5 6, valuelabel labsize(medsmall)) ylabel(,labsize(large))) ///
						(connected SHOULD_GAS_BUT_HAS_HYBRID hhfaminc_2017_groups, lpattern(-) lcolor(green) mcolor(green) msymbol(D) msize(large) lwidth(1) xlabel(1 2 3 4 5 6, valuelabel labsize(medsmall)) ylabel(,labsize(large))) ///
						(rarea hi1 lo1 hhfaminc_2017_groups, lcolor(grey%30) color(grey%30)) ///
						(rarea hi2 lo2 hhfaminc_2017_groups, lcolor(grey%30) color(grey%30)) ///
						, ///
				   graphregion(color(white)) /*title("Likely Mistakes by Demographic")*/ ///
					xtitle("Income group", size(vlarge))  ytitle("Share of group (pp. difference)", size(vlarge)) ///
					legend( label(1 "p(Gas | High Expenses)") label(2 "p(Hybrid | Low Expenses)")  ///
						 position(7) bmargin(medium) symxsize(medium) size(medium) rows(1) ring(0) region(lwidth(none)) ///
						 order(2 1 )) 	
				
				graph export "output/Figure_04b.png", as(png) width(1000) replace				
			restore
		
		
*========================================================================
**# Figure 5
*========================================================================
		
	*== DEFINING DEMOGRAPHIC VARIABLES
			*binary indicator for male
			capture drop MALE
			gen MALE = 2-r_sex_imp
			replace MALE = 2-r_sex if MALE == .
			*age
			capture drop R_AGE
			gen R_AGE=r_age
			replace R_AGE = r_age_imp if R_AGE <=0
				capture drop R_AGE_group
				gen R_AGE_group = .
				replace R_AGE_group = 1 if R_AGE < 40 & R_AGE > 0
				replace R_AGE_group = 2 if R_AGE < 60 & R_AGE >= 40
				replace R_AGE_group = 3 if R_AGE < 100 & R_AGE >= 60
			*rural area
			capture drop RURAL
			gen RURAL = 0 if urban !=.
			replace RURAL = 1 if urban == 4
			*educational attainment
			capture drop EDUC
			gen EDUC = educ
			*educational attainment
			capture drop EDUC
			gen EDUC = educ
			capture drop EDUC_group
			gen EDUC_group = .
			replace EDUC_group = 1 if EDUC == 1 | EDUC == 2
			replace EDUC_group = 2 if EDUC == 3 | EDUC == 4
			replace EDUC_group = 3 if EDUC == 5
 
	 program define demographics 

		sum SHOULD_GAS_BUT_HAS_HYBRID if MALE == 0
			scalar female_H_mean = r(mean)
			scalar female_H_sd = r(sd)
			scalar N = r(N)
			scalar female_H_sd = female_H_sd/sqrt(N)
			sum SHOULD_GAS_BUT_HAS_HYBRID if MALE == 1
			scalar male_H_mean = r(mean)
			scalar male_H_sd = r(sd)
			scalar N = r(N)
			scalar male_H_sd = male_H_sd/sqrt(N)
			sum SHOULD_HYBRID_BUT_HAS_GAS if MALE == 0
			scalar female_NH_mean = r(mean)
			scalar female_NH_sd = r(sd)
			scalar N = r(N)
			scalar female_NH_sd = female_NH_sd/sqrt(N)
			sum SHOULD_HYBRID_BUT_HAS_GAS if MALE == 1
			scalar male_NH_mean = r(mean)
			scalar male_NH_sd = r(sd)
			scalar N = r(N)
			scalar male_NH_sd = male_NH_sd/sqrt(N)

			sum SHOULD_GAS_BUT_HAS_HYBRID if R_AGE < 40 & R_AGE > 0
			scalar age0to40_H_mean = r(mean)
			scalar age0to40_H_sd = r(sd)
			scalar N = r(N)
			scalar age0to40_H_sd = age0to40_H_sd/sqrt(N)
			sum SHOULD_HYBRID_BUT_HAS_GAS if R_AGE < 40 & R_AGE > 0
			scalar age0to40_NH_mean = r(mean)
			scalar age0to40_NH_sd = r(sd)
			scalar N = r(N)
			scalar age0to40_NH_sd = age0to40_NH_sd/sqrt(N)
			sum SHOULD_GAS_BUT_HAS_HYBRID if R_AGE < 60 & R_AGE >= 40
			scalar age40to60_H_mean = r(mean)
			scalar age40to60_H_sd = r(sd)
			scalar N = r(N)
			scalar age40to60_H_sd = age40to60_H_sd/sqrt(N)
			sum SHOULD_HYBRID_BUT_HAS_GAS if R_AGE < 60 & R_AGE >= 40
			scalar age40to60_NH_mean = r(mean)
			scalar age40to60_NH_sd = r(sd)
			scalar N = r(N)
			scalar age40to60_NH_sd = age40to60_NH_sd/sqrt(N)
			sum SHOULD_GAS_BUT_HAS_HYBRID if R_AGE < 100 & R_AGE >= 60
			scalar age60to100_H_mean = r(mean)
			scalar age60to100_H_sd = r(sd)
			scalar N = r(N)
			scalar age60to100_H_sd = age60to100_H_sd/sqrt(N)
			sum SHOULD_HYBRID_BUT_HAS_GAS if R_AGE < 100 & R_AGE >= 60
			scalar age60to100_NH_mean = r(mean)
			scalar age60to100_NH_sd = r(sd)
			scalar N = r(N)
			scalar age60to100_NH_sd = age60to100_NH_sd/sqrt(N)

			sum SHOULD_GAS_BUT_HAS_HYBRID if RURAL == 0
			scalar urban_H_mean = r(mean)
			scalar urban_H_sd = r(sd)
			scalar N = r(N)
			scalar urban_H_sd = urban_H_sd/sqrt(N)
			sum SHOULD_GAS_BUT_HAS_HYBRID if RURAL == 1
			scalar rural_H_mean = r(mean)
			scalar rural_H_sd = r(sd)
			scalar N = r(N)
			scalar rural_H_sd = rural_H_sd/sqrt(N)
			sum SHOULD_HYBRID_BUT_HAS_GAS if RURAL == 0
			scalar urban_NH_mean = r(mean)
			scalar urban_NH_sd = r(sd)
			scalar N = r(N)
			scalar urban_NH_sd = urban_NH_sd/sqrt(N)
			sum SHOULD_HYBRID_BUT_HAS_GAS if RURAL == 1
			scalar rural_NH_mean = r(mean)
			scalar rural_NH_sd = r(sd)
			scalar N = r(N)
			scalar rural_NH_sd = rural_NH_sd/sqrt(N)

			sum SHOULD_GAS_BUT_HAS_HYBRID if EDUC == 1 | EDUC == 2
			scalar HSorless_H_mean = r(mean)
			scalar HSorless_H_sd = r(sd)
			scalar N = r(N)
			scalar HSorless_H_sd = HSorless_H_sd/sqrt(N)
			sum SHOULD_HYBRID_BUT_HAS_GAS if EDUC == 1 | EDUC == 2
			scalar HSorless_NH_mean = r(mean)
			scalar HSorless_NH_sd = r(sd)
			scalar N = r(N)
			scalar HSorless_NH_sd = HSorless_NH_sd/sqrt(N)
			sum SHOULD_GAS_BUT_HAS_HYBRID if EDUC == 3 | EDUC == 4
			scalar somecollege_H_mean = r(mean)
			scalar somecollege_H_sd = r(sd)
			scalar N = r(N)
			scalar somecollege_H_sd = somecollege_H_sd/sqrt(N)
			sum SHOULD_HYBRID_BUT_HAS_GAS if EDUC == 3 | EDUC == 4
			scalar somecollege_NH_mean = r(mean)
			scalar somecollege_NH_sd = r(sd)
			scalar N = r(N)
			scalar somecollege_NH_sd = somecollege_NH_sd/sqrt(N)
			sum SHOULD_GAS_BUT_HAS_HYBRID if EDUC == 5
			scalar graduate_H_mean = r(mean)
			scalar graduate_H_sd = r(sd)
			scalar N = r(N)
			scalar graduate_H_sd = graduate_H_sd/sqrt(N)
			sum SHOULD_HYBRID_BUT_HAS_GAS if EDUC == 5
			scalar graduate_NH_mean = r(mean)
			scalar graduate_NH_sd = r(sd)
			scalar N = r(N)
			scalar graduate_NH_sd = graduate_NH_sd/sqrt(N)

			matrix Hdraw = J(10,3,.)
			matrix NHdraw = J(10,3,.)
			matrix Hdraw[1,1]=male_H_mean
			matrix Hdraw[1,2]=male_H_mean - 1.96*male_H_sd
			matrix Hdraw[1,3]=male_H_mean + 1.96*male_H_sd
			matrix Hdraw[2,1]=female_H_mean
			matrix Hdraw[2,2]=female_H_mean - 1.96*female_H_sd
			matrix Hdraw[2,3]=female_H_mean + 1.96*female_H_sd
			matrix NHdraw[1,1]=male_NH_mean
			matrix NHdraw[1,2]=male_NH_mean - 1.96*male_NH_sd
			matrix NHdraw[1,3]=male_NH_mean + 1.96*male_NH_sd
			matrix NHdraw[2,1]=female_NH_mean
			matrix NHdraw[2,2]=female_NH_mean - 1.96*female_NH_sd
			matrix NHdraw[2,3]=female_NH_mean + 1.96*female_NH_sd
			matrix Hdraw[3,1]=age0to40_H_mean
			matrix Hdraw[3,2]=age0to40_H_mean - 1.96*age0to40_H_sd
			matrix Hdraw[3,3]=age0to40_H_mean + 1.96*age0to40_H_sd
			matrix NHdraw[3,1]=age0to40_NH_mean
			matrix NHdraw[3,2]=age0to40_NH_mean - 1.96*age0to40_NH_sd
			matrix NHdraw[3,3]=age0to40_NH_mean + 1.96*age0to40_NH_sd
			matrix Hdraw[4,1]=age40to60_H_mean
			matrix Hdraw[4,2]=age40to60_H_mean - 1.96*age40to60_H_sd
			matrix Hdraw[4,3]=age40to60_H_mean + 1.96*age40to60_H_sd
			matrix NHdraw[4,1]=age40to60_NH_mean
			matrix NHdraw[4,2]=age40to60_NH_mean - 1.96*age40to60_NH_sd
			matrix NHdraw[4,3]=age40to60_NH_mean + 1.96*age40to60_NH_sd
			matrix Hdraw[5,1]=age60to100_H_mean
			matrix Hdraw[5,2]=age60to100_H_mean - 1.96*age60to100_H_sd
			matrix Hdraw[5,3]=age60to100_H_mean + 1.96*age60to100_H_sd
			matrix NHdraw[5,1]=age60to100_NH_mean
			matrix NHdraw[5,2]=age60to100_NH_mean - 1.96*age60to100_NH_sd
			matrix NHdraw[5,3]=age60to100_NH_mean + 1.96*age60to100_NH_sd
			matrix Hdraw[6,1]=rural_H_mean
			matrix Hdraw[6,2]=rural_H_mean - 1.96*rural_H_sd
			matrix Hdraw[6,3]=rural_H_mean + 1.96*rural_H_sd
			matrix Hdraw[7,1]=urban_H_mean
			matrix Hdraw[7,2]=urban_H_mean - 1.96*urban_H_sd
			matrix Hdraw[7,3]=urban_H_mean + 1.96*urban_H_sd
			matrix NHdraw[6,1]=rural_NH_mean
			matrix NHdraw[6,2]=rural_NH_mean - 1.96*rural_NH_sd
			matrix NHdraw[6,3]=rural_NH_mean + 1.96*rural_NH_sd
			matrix NHdraw[7,1]=urban_NH_mean
			matrix NHdraw[7,2]=urban_NH_mean - 1.96*urban_NH_sd
			matrix NHdraw[7,3]=urban_NH_mean + 1.96*urban_NH_sd
			matrix Hdraw[8,1]=HSorless_H_mean
			matrix Hdraw[8,2]=HSorless_H_mean - 1.96*HSorless_H_sd
			matrix Hdraw[8,3]=HSorless_H_mean + 1.96*HSorless_H_sd
			matrix NHdraw[8,1]=HSorless_NH_mean
			matrix NHdraw[8,2]=HSorless_NH_mean - 1.96*HSorless_NH_sd
			matrix NHdraw[8,3]=HSorless_NH_mean + 1.96*HSorless_NH_sd
			matrix Hdraw[9,1]=somecollege_H_mean
			matrix Hdraw[9,2]=somecollege_H_mean - 1.96*somecollege_H_sd
			matrix Hdraw[9,3]=somecollege_H_mean + 1.96*somecollege_H_sd
			matrix NHdraw[9,1]=somecollege_NH_mean
			matrix NHdraw[9,2]=somecollege_NH_mean - 1.96*somecollege_NH_sd
			matrix NHdraw[9,3]=somecollege_NH_mean + 1.96*somecollege_NH_sd
			matrix Hdraw[10,1]=graduate_H_mean
			matrix Hdraw[10,2]=graduate_H_mean - 1.96*graduate_H_sd
			matrix Hdraw[10,3]=graduate_H_mean + 1.96*graduate_H_sd
			matrix NHdraw[10,1]=graduate_NH_mean
			matrix NHdraw[10,2]=graduate_NH_mean - 1.96*graduate_NH_sd
			matrix NHdraw[10,3]=graduate_NH_mean + 1.96*graduate_NH_sd
		end
			
			
		preserve		
			demographics
				matrix Hdraw_diff = Hdraw
				matrix NHdraw_diff = NHdraw

				quietly sum HAS
				scalar mean = r(mean)								
				forvalues i=1/10 {
				matrix Hdraw_diff[`i',1]=Hdraw[`i',1] - mean
				matrix Hdraw_diff[`i',2]=Hdraw[`i',2] - mean
				matrix Hdraw_diff[`i',3]=Hdraw[`i',3] - mean
				matrix NHdraw_diff[`i',1]=NHdraw[`i',1] - 1 + mean
				matrix NHdraw_diff[`i',2]=NHdraw[`i',2] - 1 + mean
				matrix NHdraw_diff[`i',3]=NHdraw[`i',3] - 1 + mean
				}	

			coefplot (matrix(Hdraw_diff[,1]), msymbol(D) ci((Hdraw_diff[,2] Hdraw_diff[,3])) label(p(Hybrid | Low Expenses))) (matrix(NHdraw_diff[,1]), msymbol(Oh) ci((NHdraw_diff[,2] NHdraw_diff[,3])) label(p(Gas | High Expenses))), scheme(s1color) rescale(100) xtitle("{bf:Share of group (pp. difference)}") groups(r1 r2 = "{bf:Gender}" r3 r4 r5= "{bf:Age}" r6 r7= "{bf:Location}" r8 r9 r10= "{bf:Education}") coeflabels(r1 = "Male" r2 = "Female"  r3= "< 40" r4= "40 to 60" r5 = "> 60"  r6 = "Rural" r7 = "Urban" r8 = "HS or less" r9 = "(Some) College" r10 = "Graduate Degree")
			graph export "output/Figure_05.png", as(png) width(1000) replace	
		restore
			


*========================================================================
**# Table 3
*========================================================================
* NOTE: These vehicle characteristics from Wards Automotive are subject to use restrictions: length width height weight liters valves horsepower rpm
* They can be merged in after gaining data access (see Readme.txt) to exactly reproduce this table.
/*		preserve
			capture drop HASNT
			capture gen HASNT = 1-HAS
			replace HASNT = . if HAS == .
			* output coefficients per 1000 USD
			scalar norm = 1/1000
			replace PricePremium_initial_real = PricePremium_initial_real * norm
			capture drop CumulativeFuelSavings CumulativeFuelSavings_3pc CumulativeFuelSavings_7pc
			gen CumulativeFuelSavings_3pc = InitialSavings_real * 11.30 * norm
			gen CumulativeFuelSavings_7pc = InitialSavings_real * 8.75 * norm
		
			capture egen vehyear_vehtype = group(vehyear vehtype)

			capture drop CumulativeFuelSavings
			gen CumulativeFuelSavings = CumulativeFuelSavings_3pc
			capture drop UpfrontInvestment
			gen UpfrontInvestment = PricePremium_initial_real
			
			* P(HYBRID|LOWCOST) [3% discounting]
				reghdfe HAS  i.hhfaminc_2017_groups i.EDUC_group i.R_AGE_group i.MALE i.urbrur length width height weight liters valves horsepower rpm if APPROACH1_LO_SHOULD ==0, absorb(vehyear_vehtype maketextNHTS hhstate)			
					outreg2 using "output/Table_03.doc" , alpha(.05) symbol(*) slow(1000) replace ctitle(EPA_3pc) addtext(State FE, YES, Year-Type FE, YES, Make FE, YES, Model FE, NO, Engineering Specs, YES, Fuel Economy, MPG_ideal)
							
				reghdfe HAS CumulativeFuelSavings UpfrontInvestment i.hhfaminc_2017_groups i.EDUC_group i.R_AGE_group i.MALE i.urbrur length width height weight liters valves horsepower rpm if APPROACH1_LO_SHOULD ==0, absorb(vehyear_vehtype maketextNHTS hhstate)			
					quietly scalar gamma_hat = - _b[CumulativeFuelSavings] / _b[UpfrontInvestment]	
					outreg2 using "output/Table_03.doc" , alpha(.05) symbol(*) slow(1000) append ctitle(EPA_3pc) addstat(GammaHat, gamma_hat) addtext(State FE, YES, Year-Type FE, YES, Make FE, YES, Model FE, NO, Engineering Specs, YES, Fuel Economy, MPG_ideal)
					capture gen over50k  = 0 if hhfaminc_2017_groups !=.
					capture replace over50k = 1 if hhfaminc_2017_groups == 3 | hhfaminc_2017_groups == 4 | hhfaminc_2017_groups == 5 | hhfaminc_2017_groups == 6
					capture gen Savings_over50k = over50k * CumulativeFuelSavings
				reghdfe HAS CumulativeFuelSavings Savings_over50k UpfrontInvestment i.hhfaminc_2017_groups i.EDUC_group i.R_AGE_group i.MALE i.urbrur length width height weight liters valves horsepower rpm  if APPROACH1_LO_SHOULD ==0, absorb(vehyear_vehtype maketextNHTS hhstate)			
				quietly scalar gamma_hat = - _b[CumulativeFuelSavings] / _b[UpfrontInvestment]	
					quietly scalar gamma_hat50k = gamma_hat - _b[Savings_over50k] / _b[UpfrontInvestment]	
							outreg2 using "output/Table_03.doc" , alpha(.05) symbol(*) slow(1000) append ctitle(EPA_3pc)  addstat(GammaHat (under 50k), gamma_hat, GammaHat (over 50k), gamma_hat50k) addtext(State FE, YES, Year-Type FE, YES, Make FE, YES, Model FE, NO, Engineering Specs, YES, Fuel Economy, MPG_ideal)
		
			* P(GAS|LOWCOST) [3% discounting]		
				reghdfe HASNT i.hhfaminc_2017_groups i.EDUC_group i.R_AGE_group i.MALE i.urbrur length width height weight liters valves horsepower rpm if APPROACH1_LO_SHOULD ==1, absorb(vehyear_vehtype maketextNHTS hhstate)			
					outreg2 using "output/Table_03.doc" , alpha(.05) symbol(*) slow(1000) append ctitle(EPA_3pc)  addtext(State FE, YES, Year-Type FE, YES, Make FE, YES, Model FE, NO, Engineering Specs, YES, Fuel Economy, MPG_ideal)
				reghdfe HASNT CumulativeFuelSavings UpfrontInvestment i.hhfaminc_2017_groups i.EDUC_group i.R_AGE_group i.MALE i.urbrur length width height weight liters valves horsepower rpm if APPROACH1_LO_SHOULD ==1, absorb(vehyear_vehtype maketextNHTS hhstate)			
					quietly scalar gamma_hat = - _b[CumulativeFuelSavings] / _b[UpfrontInvestment]	
					outreg2 using "output/Table_03.doc" , alpha(.05) symbol(*) slow(1000) append ctitle(EPA_3pc)  addstat(GammaHat, gamma_hat) addtext(State FE, YES, Year-Type FE, YES, Make FE, YES, Model FE, NO, Engineering Specs, YES, Fuel Economy, MPG_ideal)
					capture gen over50k  = 0 if hhfaminc_2017_groups !=.
					capture replace over50k = 1 if hhfaminc_2017_groups == 3 | hhfaminc_2017_groups == 4 | hhfaminc_2017_groups == 5 | hhfaminc_2017_groups == 6
					capture gen Savings_over50k = over50k * CumulativeFuelSavings
				reghdfe HASNT CumulativeFuelSavings Savings_over50k UpfrontInvestment i.hhfaminc_2017_groups i.EDUC_group i.R_AGE_group i.MALE i.urbrur length width height weight liters valves horsepower rpm  if APPROACH1_LO_SHOULD ==1, absorb(vehyear_vehtype maketextNHTS hhstate)			
					quietly scalar gamma_hat = - _b[CumulativeFuelSavings] / _b[UpfrontInvestment]	
					quietly scalar gamma_hat50k = gamma_hat - _b[Savings_over50k] / _b[UpfrontInvestment]	
							outreg2 using "output/Table_03.doc" , alpha(.05) symbol(*) slow(1000) append ctitle(EPA_3pc)  addstat(GammaHat (under 50k), gamma_hat, GammaHat (over 50k), gamma_hat50k) addtext(State FE, YES, Year-Type FE, YES, Make FE, YES, Model FE, NO, Engineering Specs, YES, Fuel Economy, MPG_ideal)
		
		restore
*/
			
*========================================================================
**# Table 6
*========================================================================
		* Panel (a): Same approach as in Table 2, but relies on restricted access data from MaritzCX
		* Panel (b): Same approach as in Table 2, but using TrueCar used car price data
			gen APPROACH5_HI_SHOULD = .
			replace APPROACH5_HI_SHOULD = 0
			replace APPROACH5_HI_SHOULD = 1 if FUELSAVING_vs_PREMIUM_used*(14-vehage)/14 > APPROACH1_HI_cutoff
			replace APPROACH5_HI_SHOULD = . if FUELSAVING_vs_PREMIUM_used ==.
			gen APPROACH5_LO_SHOULD = .
			replace APPROACH5_LO_SHOULD = 0
			replace APPROACH5_LO_SHOULD = 1 if FUELSAVING_vs_PREMIUM_used*(14-vehage)/14 > APPROACH1_LO_cutoff
			replace APPROACH5_LO_SHOULD = . if FUELSAVING_vs_PREMIUM_used ==.
	
		preserve
		keep if APPROACH5_LO_SHOULD !=.
		outreg2 HAS APPROACH5_LO_SHOULD using "output/Table_06b.xls", append cross  noaster
		restore			
		preserve
		keep if APPROACH5_LO_SHOULD !=.
		outreg2 HAS APPROACH5_HI_SHOULD using "output/Table_06b.xls", append cross  noaster
		restore		
			
				
 
*========================================================================
* SWITCH TO ALL CAR SAMPLE
*========================================================================

use "VEH2.dta", clear
		* STRICT CUTOFF DISCOUNT RATES: REPLICATE EPA CAFE 2017-2015 ASSESSMENT (3% and 7% over 14 years)
		* cutoffs are calculated for an average 14 years of discounted fuel savings 
		scalar APPROACH1_HI_cutoff = 1/8.75
		scalar APPROACH1_LO_cutoff = 1/11.30	
		capture drop MPGsavings_EPA_3pc MPGsavings_EPA_7pc
		gen MPGsavings_EPA_3pc = MPGcost / APPROACH1_LO_cutoff
		gen MPGsavings_EPA_7pc = MPGcost / APPROACH1_HI_cutoff		
		scalar MPGprice = 106
		scalar MPGprice_HI = 207

	* IDENTIFY DRIVERS WITH "TOO LITTLE" MPG
		gen EPA_3pc_toolittle = 0 if MPGsavings_EPA_3pc !=.
		replace EPA_3pc_toolittle = 1 if MPGsavings_EPA_3pc >= MPGprice & MPGsavings_EPA_3pc !=.
		gen EPA_7pc_toolittle = 0 if MPGsavings_EPA_7pc !=.
		replace EPA_7pc_toolittle = 1 if MPGsavings_EPA_7pc >= MPGprice_HI & MPGsavings_EPA_7pc !=.
		* OPTIMAL MPG: mu* = (VMT * Pgas)/sqrt(Pmpg * F)
		capture drop mpg_ideal*
		gen mpg_ideal_3pc  = ((bestmile * EIA_gasprice_surveyyear_real) / (MPGprice * APPROACH1_LO_cutoff))^0.5
		gen mpg_ideal_7pc  = ((bestmile * EIA_gasprice_surveyyear_real) / (MPGprice_HI * APPROACH1_HI_cutoff))^0.5
		* again, but using median MPG by NHTS wave
		capture drop RND_mpg_Wards_combined
		gen RND_mpg_Wards_combined = NHTSwave_mpgmedian
		gen RND_MPGcost = (bestmile / (RND_mpg_Wards_combined)^2) * EIA_gasprice_surveyyear_real	
		gen RND_MPGsavings_EPA_3pc = RND_MPGcost / APPROACH1_LO_cutoff
		gen RND_MPGsavings_EPA_7pc = RND_MPGcost / APPROACH1_HI_cutoff		
		gen RND_EPA_3pc_toolittle = 0 if RND_MPGsavings_EPA_3pc !=.
		replace RND_EPA_3pc_toolittle = 1 if RND_MPGsavings_EPA_3pc >= MPGprice & RND_MPGsavings_EPA_3pc !=.
		gen RND_EPA_7pc_toolittle = 0 if RND_MPGsavings_EPA_7pc !=.
		replace RND_EPA_7pc_toolittle = 1 if RND_MPGsavings_EPA_7pc >= MPGprice_HI & RND_MPGsavings_EPA_7pc !=.
		
	
*========================================================================
**# Table 4
*========================================================================
		preserve
			keep if EPA_3pc_toolittle !=.
			outreg2 NHTSwave EPA_3pc_toolittle  using  "output/Table_04.xls", replace cross noaster
		restore
		preserve
			keep if EPA_7pc_toolittle !=.
			outreg2 NHTSwave EPA_7pc_toolittle  using  "output/Table_04.xls", append cross noaster
		restore		
				
				/* mu* and foregone savings were pre-calculated before dropping mpg_wards for Dataverse using the following code:
				capture drop mpg_ideal*
				gen mpg_ideal_3pc  = ((bestmile * EIA_gasprice_surveyyear_real) / (MPGprice * APPROACH1_LO_cutoff))^0.5
				gen mpg_ideal_7pc  = ((bestmile * EIA_gasprice_surveyyear_real) / (MPGprice_HI * APPROACH1_HI_cutoff))^0.5
				capture drop FOREGONE_*
				gen FOREGONE_fuelsavings_3pc = bestmile * EIA_gasprice_surveyyear_real * ((1/mpg_Wards_combined) - (1/mpg_ideal_3pc)) / APPROACH1_LO_cutoff
				gen FOREGONE_vehiclecost_3pc = MPGprice * (mpg_ideal_3pc - mpg_Wards_combined)
				gen FOREGONE_savings_3pc = FOREGONE_fuelsavings_3pc - FOREGONE_vehiclecost_3pc
				gen FOREGONE_fuelsavings_7pc = bestmile * EIA_gasprice_surveyyear_real * ((1/mpg_Wards_combined) - (1/mpg_ideal_7pc)) / APPROACH1_HI_cutoff
				gen FOREGONE_vehiclecost_7pc = MPGprice_HI * (mpg_ideal_7pc - mpg_Wards_combined)
				gen FOREGONE_savings_7pc = FOREGONE_fuelsavings_7pc - FOREGONE_vehiclecost_7pc
				*/			  				  
				sum FOREGONE_savings_3pc if EPA_3pc_toolittle == 0
				sum FOREGONE_savings_3pc if EPA_3pc_toolittle == 1
				sum FOREGONE_savings_7pc if EPA_7pc_toolittle == 0 
				sum FOREGONE_savings_7pc if EPA_7pc_toolittle == 1
		
				
*========================================================================
**# Figure 6
*========================================================================
		capture drop MPGtoolow RND_MPGtoolow
		gen MPGtoolow = EPA_3pc_toolittle
		gen RND_MPGtoolow = RND_EPA_3pc_toolittle							

		*-----------------------------------------------------------------------
		* Panel (a): absolute shares
		*-----------------------------------------------------------------------
		preserve
		keep if NHTSwave==2017
		replace MPGtoolow = MPGtoolow * 100
		collapse (mean) MPGtoolow (semean) stdev=MPGtoolow, by(hhfaminc_2017_groups)
			gen hi1=MPGtoolow+1.96*stdev
			gen lo1=MPGtoolow-1.96*stdev
			twoway (connected MPGtoolow hhfaminc_2017_groups, lwidth(1) lcolor(red) mcolor(red) msymbol(Oh) msize(large) lwidth(1) xlabel(1 2 3 4 5 6, valuelabel labsize(medsmall)) ylabel(,labsize(large)))  ///
					(rarea hi1 lo1 hhfaminc_2017_groups, lcolor(grey%30) color(grey%30)) ///
					, ///
			   graphregion(color(white)) /*title("Likely Mistakes by Demographic")*/ ///
				xtitle("Income group", size(vlarge))  ytitle("Share of group (in %)", size(vlarge)) /* yscale(range(35 60)) */ ///
				legend(off) 	
 		graph export "output/Figure_06a.png", as(png) width(1000) replace				
		restore
	
 		*-----------------------------------------------------------------------
		* Panel (b): absolute share minus  relative to counterfactual where everyon has median(mpg)
		*-----------------------------------------------------------------------
		preserve
		keep if NHTSwave==2017
		replace MPGtoolow = (MPGtoolow - RND_MPGtoolow) * 100
		collapse (mean) MPGtoolow (semean) stdev=MPGtoolow, by(hhfaminc_2017_groups)
			gen hi1=MPGtoolow+1.96*stdev
			gen lo1=MPGtoolow-1.96*stdev
			twoway (connected MPGtoolow hhfaminc_2017_groups, lwidth(1) lcolor(red) mcolor(red) msymbol(Oh) msize(large) lwidth(1) xlabel(1 2 3 4 5 6, valuelabel labsize(medsmall)) ylabel(,labsize(large))) ///
					(rarea hi1 lo1 hhfaminc_2017_groups, lcolor(grey%30) color(grey%30)) ///
					, ///
			   graphregion(color(white)) /*title("Likely Mistakes by Demographic")*/ ///
				xtitle("Income group", size(vlarge))  ytitle("Share of group (pp. difference)", size(vlarge)) ///
				legend(off) 	
 		graph export "output/Figure_06b.png", as(png) width(1000) replace				
		restore
		
		
		
*========================================================================
**# Figure 7
*========================================================================
		sum MPGtoolow if MALE == 0
		scalar female_H_mean = r(mean)
		scalar female_H_sd = r(sd)
		scalar N = r(N)
		scalar female_H_sd = female_H_sd/sqrt(N)
		sum MPGtoolow if MALE == 1
		scalar male_H_mean = r(mean)
		scalar male_H_sd = r(sd)
		scalar N = r(N)
		scalar male_H_sd = male_H_sd/sqrt(N)
		sum MPGtoolow if R_AGE < 40 & R_AGE > 0
		scalar age0to40_H_mean = r(mean)
		scalar age0to40_H_sd = r(sd)
		scalar N = r(N)
		scalar age0to40_H_sd = age0to40_H_sd/sqrt(N)
		sum MPGtoolow if R_AGE < 40 & R_AGE > 0
		scalar age0to40_NH_mean = r(mean)
		scalar age0to40_NH_sd = r(sd)
		scalar N = r(N)
		scalar age0to40_NH_sd = age0to40_NH_sd/sqrt(N)
		sum MPGtoolow if R_AGE < 60 & R_AGE >= 40
		scalar age40to60_H_mean = r(mean)
		scalar age40to60_H_sd = r(sd)
		scalar N = r(N)
		scalar age40to60_H_sd = age40to60_H_sd/sqrt(N)
		sum MPGtoolow if R_AGE < 100 & R_AGE >= 60
		scalar age60to100_H_mean = r(mean)
		scalar age60to100_H_sd = r(sd)
		scalar N = r(N)
		scalar age60to100_H_sd = age60to100_H_sd/sqrt(N)
		sum MPGtoolow if RURAL == 0
		scalar urban_H_mean = r(mean)
		scalar urban_H_sd = r(sd)
		scalar N = r(N)
		scalar urban_H_sd = urban_H_sd/sqrt(N)
		sum MPGtoolow if RURAL == 1
		scalar rural_H_mean = r(mean)
		scalar rural_H_sd = r(sd)
		scalar N = r(N)
		scalar rural_H_sd = rural_H_sd/sqrt(N)
		sum MPGtoolow if EDUC == 1 | EDUC == 2
		scalar HSorless_H_mean = r(mean)
		scalar HSorless_H_sd = r(sd)
		scalar N = r(N)
		scalar HSorless_H_sd = HSorless_H_sd/sqrt(N)
		sum MPGtoolow if EDUC == 3 | EDUC == 4
		scalar somecollege_H_mean = r(mean)
		scalar somecollege_H_sd = r(sd)
		scalar N = r(N)
		scalar somecollege_H_sd = somecollege_H_sd/sqrt(N)
		sum MPGtoolow if EDUC == 5
		scalar graduate_H_mean = r(mean)
		scalar graduate_H_sd = r(sd)
		scalar N = r(N)
		scalar graduate_H_sd = graduate_H_sd/sqrt(N)
		matrix Hdraw = J(10,3,.)
		matrix Hdraw[1,1]=male_H_mean
		matrix Hdraw[1,2]=male_H_mean - 1.96*male_H_sd
		matrix Hdraw[1,3]=male_H_mean + 1.96*male_H_sd
		matrix Hdraw[2,1]=female_H_mean
		matrix Hdraw[2,2]=female_H_mean - 1.96*female_H_sd
		matrix Hdraw[2,3]=female_H_mean + 1.96*female_H_sd
		matrix Hdraw[3,1]=age0to40_H_mean
		matrix Hdraw[3,2]=age0to40_H_mean - 1.96*age0to40_H_sd
		matrix Hdraw[3,3]=age0to40_H_mean + 1.96*age0to40_H_sd
		matrix Hdraw[4,1]=age40to60_H_mean
		matrix Hdraw[4,2]=age40to60_H_mean - 1.96*age40to60_H_sd
		matrix Hdraw[4,3]=age40to60_H_mean + 1.96*age40to60_H_sd
		matrix Hdraw[5,1]=age60to100_H_mean
		matrix Hdraw[5,2]=age60to100_H_mean - 1.96*age60to100_H_sd
		matrix Hdraw[5,3]=age60to100_H_mean + 1.96*age60to100_H_sd
		matrix Hdraw[6,1]=rural_H_mean
		matrix Hdraw[6,2]=rural_H_mean - 1.96*rural_H_sd
		matrix Hdraw[6,3]=rural_H_mean + 1.96*rural_H_sd
		matrix Hdraw[7,1]=urban_H_mean
		matrix Hdraw[7,2]=urban_H_mean - 1.96*urban_H_sd
		matrix Hdraw[7,3]=urban_H_mean + 1.96*urban_H_sd
		matrix Hdraw[8,1]=HSorless_H_mean
		matrix Hdraw[8,2]=HSorless_H_mean - 1.96*HSorless_H_sd
		matrix Hdraw[8,3]=HSorless_H_mean + 1.96*HSorless_H_sd
		matrix Hdraw[9,1]=somecollege_H_mean
		matrix Hdraw[9,2]=somecollege_H_mean - 1.96*somecollege_H_sd
		matrix Hdraw[9,3]=somecollege_H_mean + 1.96*somecollege_H_sd
		matrix Hdraw[10,1]=graduate_H_mean
		matrix Hdraw[10,2]=graduate_H_mean - 1.96*graduate_H_sd
		matrix Hdraw[10,3]=graduate_H_mean + 1.96*graduate_H_sd

				sum RND_MPGtoolow if MALE == 0
				scalar female_H_mean = r(mean)
				scalar female_H_sd = r(sd)
				scalar N = r(N)
				scalar female_H_sd = female_H_sd/sqrt(N)
				sum RND_MPGtoolow if MALE == 1
				scalar male_H_mean = r(mean)
				scalar male_H_sd = r(sd)
				scalar N = r(N)
				scalar male_H_sd = male_H_sd/sqrt(N)
				sum RND_MPGtoolow if R_AGE < 40 & R_AGE > 0
				scalar age0to40_H_mean = r(mean)
				scalar age0to40_H_sd = r(sd)
				scalar N = r(N)
				scalar age0to40_H_sd = age0to40_H_sd/sqrt(N)
				sum RND_MPGtoolow if R_AGE < 40 & R_AGE > 0
				scalar age0to40_NH_mean = r(mean)
				scalar age0to40_NH_sd = r(sd)
				scalar N = r(N)
				scalar age0to40_NH_sd = age0to40_NH_sd/sqrt(N)
				sum RND_MPGtoolow if R_AGE < 60 & R_AGE >= 40
				scalar age40to60_H_mean = r(mean)
				scalar age40to60_H_sd = r(sd)
				scalar N = r(N)
				scalar age40to60_H_sd = age40to60_H_sd/sqrt(N)
				sum RND_MPGtoolow if R_AGE < 100 & R_AGE >= 60
				scalar age60to100_H_mean = r(mean)
				scalar age60to100_H_sd = r(sd)
				scalar N = r(N)
				scalar age60to100_H_sd = age60to100_H_sd/sqrt(N)
				sum RND_MPGtoolow if RURAL == 0
				scalar urban_H_mean = r(mean)
				scalar urban_H_sd = r(sd)
				scalar N = r(N)
				scalar urban_H_sd = urban_H_sd/sqrt(N)
				sum RND_MPGtoolow if RURAL == 1
				scalar rural_H_mean = r(mean)
				scalar rural_H_sd = r(sd)
				scalar N = r(N)
				scalar rural_H_sd = rural_H_sd/sqrt(N)
				sum RND_MPGtoolow if EDUC == 1 | EDUC == 2
				scalar HSorless_H_mean = r(mean)
				scalar HSorless_H_sd = r(sd)
				scalar N = r(N)
				scalar HSorless_H_sd = HSorless_H_sd/sqrt(N)
				sum RND_MPGtoolow if EDUC == 3 | EDUC == 4
				scalar somecollege_H_mean = r(mean)
				scalar somecollege_H_sd = r(sd)
				scalar N = r(N)
				scalar somecollege_H_sd = somecollege_H_sd/sqrt(N)
				sum RND_MPGtoolow if EDUC == 5
				scalar graduate_H_mean = r(mean)
				scalar graduate_H_sd = r(sd)
				scalar N = r(N)
				scalar graduate_H_sd = graduate_H_sd/sqrt(N)
				matrix RND_Hdraw = J(10,3,.)
				matrix RND_Hdraw[1,1]=male_H_mean
				matrix RND_Hdraw[1,2]=male_H_mean - 1.96*male_H_sd
				matrix RND_Hdraw[1,3]=male_H_mean + 1.96*male_H_sd
				matrix RND_Hdraw[2,1]=female_H_mean
				matrix RND_Hdraw[2,2]=female_H_mean - 1.96*female_H_sd
				matrix RND_Hdraw[2,3]=female_H_mean + 1.96*female_H_sd
				matrix RND_Hdraw[3,1]=age0to40_H_mean
				matrix RND_Hdraw[3,2]=age0to40_H_mean - 1.96*age0to40_H_sd
				matrix RND_Hdraw[3,3]=age0to40_H_mean + 1.96*age0to40_H_sd
				matrix RND_Hdraw[4,1]=age40to60_H_mean
				matrix RND_Hdraw[4,2]=age40to60_H_mean - 1.96*age40to60_H_sd
				matrix RND_Hdraw[4,3]=age40to60_H_mean + 1.96*age40to60_H_sd
				matrix RND_Hdraw[5,1]=age60to100_H_mean
				matrix RND_Hdraw[5,2]=age60to100_H_mean - 1.96*age60to100_H_sd
				matrix RND_Hdraw[5,3]=age60to100_H_mean + 1.96*age60to100_H_sd
				matrix RND_Hdraw[6,1]=rural_H_mean
				matrix RND_Hdraw[6,2]=rural_H_mean - 1.96*rural_H_sd
				matrix RND_Hdraw[6,3]=rural_H_mean + 1.96*rural_H_sd
				matrix RND_Hdraw[7,1]=urban_H_mean
				matrix RND_Hdraw[7,2]=urban_H_mean - 1.96*urban_H_sd
				matrix RND_Hdraw[7,3]=urban_H_mean + 1.96*urban_H_sd
				matrix RND_Hdraw[8,1]=HSorless_H_mean
				matrix RND_Hdraw[8,2]=HSorless_H_mean - 1.96*HSorless_H_sd
				matrix RND_Hdraw[8,3]=HSorless_H_mean + 1.96*HSorless_H_sd
				matrix RND_Hdraw[9,1]=somecollege_H_mean
				matrix RND_Hdraw[9,2]=somecollege_H_mean - 1.96*somecollege_H_sd
				matrix RND_Hdraw[9,3]=somecollege_H_mean + 1.96*somecollege_H_sd
				matrix RND_Hdraw[10,1]=graduate_H_mean
				matrix RND_Hdraw[10,2]=graduate_H_mean - 1.96*graduate_H_sd
				matrix RND_Hdraw[10,3]=graduate_H_mean + 1.96*graduate_H_sd
		
			matrix Hdraw_diff = J(10,3,.)
				forvalues i=1/10 {
				matrix Hdraw_diff[`i',1]=Hdraw[`i',1] - RND_Hdraw[`i',1]
				matrix Hdraw_diff[`i',2]=Hdraw[`i',2] - RND_Hdraw[`i',1]
				matrix Hdraw_diff[`i',3]=Hdraw[`i',3] - RND_Hdraw[`i',1]
				}	
					
		coefplot (matrix(Hdraw_diff[,1]), msymbol(Oh) mcolor(red) ci((Hdraw_diff[,2] Hdraw_diff[,3])) ciopts(lcolor(red)) label(Share of drivers with too little fuel economy)), scheme(s1color) rescale(100) xtitle("{bf:Share of group with too little fuel economy (pp. difference)}") groups(r1 r2 = "{bf:Gender}" r3 r4 r5= "{bf:Age}" r6 r7= "{bf:Location}" r8 r9 r10= "{bf:Education}") coeflabels(r1 = "Male" r2 = "Female"  r3= "< 40" r4= "40 to 60" r5 = "> 60"  r6 = "Rural" r7 = "Urban" r8 = "HS or less" r9 = "(Some) College" r10 = "Graduate Degree")
 		graph export "output/Figure_07.png", as(png) width(1000) replace	
	



*========================================================================
**# Table 5
*========================================================================
* NOTE: These vehicle characteristics from Wards Automotive are subject to use restrictions: length width height weight liters valves horsepower rpm
* They can be merged in after gaining data access (see Readme.txt) to exactly reproduce this table.
/*	preserve	
			capture drop MPGsavings
			gen MPY = bestmile / 1000			
			capture egen bestmile_decile = cut(bestmile), group(10)	
			capture gen bestmile_sq = bestmile^2
			capture gen EIA_gasprice_vehicleyear_sq = EIA_gasprice_vehicleyear^2
			capture gen bestmile_gasprice = bestmile * EIA_gasprice_vehicleyear
			capture gen bestmile_gasprice_sq = bestmile_gasprice^2	
			capture egen vehyear_vehtype = group(vehyear vehtype)
									
			reghdfe EPA_3pc_toolittle i.hhfaminc_2017_groups i.EDUC_group i.R_AGE_group i.MALE i.urbrur length width height weight liters valves horsepower rpm , absorb(vehyear_vehtype maketextNHTS hhstate)			
						outreg2 using "output/Table_05.doc" , alpha(.05) symbol(*) slow(1000) dec(4) replace addtext(State FE, YES, Year-Type FE, YES, Make FE, YES, Model FE, NO, Engineering Specs, YES, Fuel Economy, MPG_ideal)
			reghdfe EPA_3pc_toolittle MPY i.hhfaminc_2017_groups i.EDUC_group i.R_AGE_group i.MALE i.urbrur length width height weight liters valves horsepower rpm , absorb(vehyear_vehtype maketextNHTS hhstate)			
						outreg2 using "output/Table_05.doc" , alpha(.05) symbol(*) slow(1000) dec(4) append addtext(State FE, YES, Year-Type FE, YES, Make FE, YES, Model FE, NO, Engineering Specs, YES, Fuel Economy, MPG_ideal)
				capture gen over50k  = 0 if hhfaminc_2017_groups !=.
				capture replace over50k = 1 if hhfaminc_2017_groups == 3 | hhfaminc_2017_groups == 4 | hhfaminc_2017_groups == 5 | hhfaminc_2017_groups == 6
				capture gen MPY_over50k = over50k * MPY
			reghdfe EPA_3pc_toolittle MPY MPY_over50k i.hhfaminc_2017_groups i.EDUC_group i.R_AGE_group i.MALE i.urbrur length width height weight liters valves horsepower rpm , absorb(vehyear_vehtype maketextNHTS hhstate)			
						outreg2 using "output/Table_05.doc" , alpha(.05) symbol(*) slow(1000) dec(4) append addtext(State FE, YES, Year-Type FE, YES, Make FE, YES, Model FE, NO, Engineering Specs, YES, Fuel Economy, MPG_ideal)				
	restore
*/	




*========================================================================
**# Table 1
*========================================================================
use "VEH1.dta", clear
*Generate the tables of means
	mean bestmile  Educ_* RURAL if HAS==0
	outreg2 bestmile Educ_* RURAL  ///
		using "output/Table_01-a.doc", replace  noaster
	capture drop tmp_*
	
	mean bestmile   Educ_* RURAL if HAS==1
	outreg2 bestmile Educ_* RURAL  ///
		using "output/Table_01-a.doc", append  noaster	
		
*Note, do income separately, since it is only for 2017 NHTS	
	mean Inc_* if HAS==0
	outreg2 Inc_*  ///
		using "output/Table_01-b.doc", replace  noaster		
	mean Inc_* if HAS==1
	outreg2 Inc_* ///
		using "output/Table_01-b.doc", append  noaster	
*Note, do age separately since it is missing in some 
	mean age if HAS==0
	outreg2 HAS age  ///
		using "output/Table_01-c.doc", replace  noaster		
	mean age if HAS==1
	outreg2 HAS age ///
		using "output/Table_01-c.doc", append  noaster	
/*Vehicle Specs (not contained in replication archive due to data restrictions)
	mean  mpg_Wards_combined length width height weight liters valves horsepower rpm vehage if HAS==0
	outreg2  mpg_Wards_combined length width height weight liters valves horsepower rpm vehage ///
		using "output/Table_01-d.doc", replace  noaster
	mean  mpg_Wards_combined  length width height weight liters valves horsepower rpm vehage if HAS==1
	outreg2 mpg_Wards_combined length width height weight liters valves horsepower rpm vehage ///
		using "output/Table_01-d.doc", append  noaster	
*/

use "VEH2.dta", clear
	mean bestmile  Educ_*  RURAL 
	outreg2 bestmile  Educ_*  RURAL ///
		using "output/Table_01-a.doc", append noaster
*Income
	mean Inc_* 
	outreg2 Inc_* using "output/Table_01-b.doc", append noaster
*Age		
	mean age 
	outreg2 age using "output/Table_01-c.doc", append noaster
/*Vehicle Specs (not contained in replication archive due to data restrictions)
	mean  mpg_Wards_combined length width height weight liters valves horsepower rpm vehage 
	outreg2  mpg_Wards_combined length width height weight liters valves horsepower rpm vehage ///
		using "output/Table_01-d.doc", append noaster
*/
